;********************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;load library for  wetbulb_stull
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/heat_stress.ncl"
;********************************************
 begin

;set baseline
 dirPath = "~/source_data/"
 dirPath2 = dirPath+"/outdoor/"

 mo  = (/"ACCESS-CM2","AWI-CM-1-1-MR","BCC-CSM2-MR","CMCC-CM2-SR5","EC-Earth3","GFDL-ESM4", \
             "IITM-ESM","KIOST-ESM","MIROC6","MIROC-ES2L","MPI-ESM1-2-HR","MRI-ESM2-0"/); 12 models 
 case = (/"ssp585" /)
 ncase = dimsizes(case);
 nmod = dimsizes(mo)

 ;start and final year of data
 syr = 1990; including extra histroical data
 fyr = 2099;
 csyr = sprinti("%0.4i",syr);
 cfyr = sprinti("%0.4i",fyr);
 nyr = fyr - syr + 1;
 ndays = nyr*365

;read dimension for 1 deg data
  dimPath=dirPath+"/domain/"
  landfrac = "sftlf_fx_360x180_1deg_gn.nc"
  area = "areacella_fx_360x180_1deg_gn.nc"

  fdim  = addfile(dimPath+landfrac, "r")
  fdim2 = addfile(dimPath+area, "r")
  lat := fdim->lat
  lon := fdim->lon
  nlat = dimsizes(lat);
  nlon = dimsizes(lon);
  garea = fdim2->areacella  ;m^2
  lfrac = fdim->sftlf/100; %
  wgt = garea*1e-6*lfrac ; sqrkm
  wgt@_FillValue = lfrac@_FillValue; make sure FillValue is set correctly for wgt
  wgt = where(wgt .eq. 0, wgt@_FillValue, wgt)
  copy_VarCoords(garea, wgt)
  printVarSummary(wgt)
  print("global total land area: "+sum(wgt)+" sq km")


  ;--*read population distribution data from "Fig.S2.Population_map_water_rural.ncl"
  ;1) population spending >30mins to access water (rescaled to UNICEF,WHO gross values)
  ;2) rural farmer population data made from 1km SMOD and 1km POP data from EU GHSL
  ;as described in Supplementary Method 2 "Population distribution data"
  dirNc =  dirPath
  filNc2  = "population_data_watercollector_farmer.nc"
  pthNc2  = dirNc+filNc2
  ncdf2 = addfile(pthNc2 ,"r")
  POPR_1deg = ncdf2->POPR_1deg; population engaged in farming
  POPR_dens = ncdf2->POPR_dens
  POPW_1deg = ncdf2->POPW_1deg; population engaged in water collection
  POPW_dens = ncdf2->POPW_dens


;===Use precalculated global warming amount relative to 1980-2009 (based on ERA5, HadCRUT5, BerkeleyEarth, NOAAGlobalTemp) to classify warming levels for each CMIP6 model
  ;see Data_global-warming-amount-ERA5-30yr-CMIP6_12models.ncl
  filename1 = dirPath+"Global_warming_amount_10yrma_"+syr+fyr+"_13models_obs30yr.csv"; 10-year running mean
  filename2 = dirPath+"Global_warming_amount_30yrma_"+syr+fyr+"_13models_obs30yr.csv"; 30-year running mean

;---Read in file as array of strings so we can parse each line
  lines  = asciiread(filename1,-1,"string")
  lines2 = asciiread(filename2,-1,"string")
  k = 1 ;number of header lines
  nlines = dimsizes(lines)-k   ; skip header lines
  delim = ","
;---Read fields (columns)
  models  =          str_get_field(lines(k:),1,delim)
  cases   =          str_get_field(lines(k:),2,delim)
  central_year10 =    str_get_field(lines(k:),3,delim)
  period10yr =       str_get_field(lines(k:),5,delim)
 ;read mean temperature anomaly from file 1
  Ta_anomaly10yr = tofloat(str_get_field(lines(k:),7,delim))
  Ta_anomaly10yr@_FillValue=1e+20

  central_year30 =    str_get_field(lines2(k:),3,delim)
  period30yr =       str_get_field(lines2(k:),5,delim)
  Ta_anomaly30yr = tofloat(str_get_field(lines2(k:),7,delim))
  Ta_anomaly30yr@_FillValue=1e+20

  ;----select either 10-yr or 30-yr period 
  Ta_anomaly = Ta_anomaly10yr;
  periods = period10yr
  central_year = central_year10


 ;==set 8 warming levels for WoE analysis but showing 7 levels in the map (or 8 with the last level ">4")
  levs = (/1,1.5,2,2.5,3,3.5,4/); 7 warming levels for map display
  nlev = dimsizes(levs)
  print("---estimate ensemble median WoE based on ensemble median Gday---")


   ;--note the ensemble median WoE map uses output from "Data.Outdoor_Gday_nomask.ncl" 

   infile = dirPath2+"Gday-gridcell-warminglevel-10yravg-ens-median-vcbc-outdoor.nc"
   fi = addfile(infile, "r")

   Gday_sun_med  = fi->Gday_sun_med
   Gday_diff_med = fi->Gday_diff_med
   Gday_zero_med = fi->Gday_zero_med
   TWday_med = fi->TWday_med
   levels = fi->lev
   printVarSummary(Gday_sun_med); [nlev,nlat,nlon]
   print(levels); 0.8 to 4.5 by 0.1 degC

   ;==find warming of emergence (WoE) of first warming amount when Gday>0 occurs at each grid cell==
   ;the 10-yr/30-yr Gday values are indexed to its middle year (e.g., 2019 for 2015-2024) corresponding to each warming amount

   WoE_Gday_sun = new((/nlat,nlon/),typeof(Gday_sun_med))
   WoE_Gday_diff = new((/nlat,nlon/),typeof(Gday_diff_med))
   WoE_Gday_zero = new((/nlat,nlon/),typeof(Gday_zero_med))
   WoE_TWday = new((/nlat,nlon/),typeof(TWday_med))
   levels := fspan(0.8,4,33); 0.8 to 4.0 by 0.1 degC
   nlev = dimsizes(levels); set WoEs>4 to NA and show by the last color bar
   do ii=0,nlev-1
     WoE_Gday_sun = where(ismissing(WoE_Gday_sun) .and. Gday_sun_med(ii,:,:) .ge. 0, levels(ii), WoE_Gday_sun)
     WoE_Gday_diff = where(ismissing(WoE_Gday_diff) .and. Gday_diff_med(ii,:,:) .ge. 0, levels(ii), WoE_Gday_diff)
     WoE_Gday_zero = where(ismissing(WoE_Gday_zero) .and. Gday_zero_med(ii,:,:) .ge. 0, levels(ii), WoE_Gday_zero)
     WoE_TWday = where(ismissing(WoE_TWday) .and. TWday_med(ii,:,:) .ge. 35, levels(ii), WoE_TWday)
   end do
   print("range of WoE_Gday_sun:")
   printMinMax(WoE_Gday_sun, 0)
   print("range of WoE_Gday_diff:")
   printMinMax(WoE_Gday_diff, 0)
   print("range of WoE_Gday_zero:")
   printMinMax(WoE_Gday_zero, 0)
   print("range of WoE_TWday:")
   printMinMax(WoE_TWday, 0)


   copy_VarCoords(Gday_sun_med(0,:,:), WoE_Gday_sun)
   copy_VarCoords(Gday_sun_med(0,:,:), WoE_Gday_diff)
   copy_VarCoords(Gday_sun_med(0,:,:), WoE_Gday_zero)
   copy_VarCoords(Gday_sun_med(0,:,:), WoE_TWday)
   printVarSummary(WoE_Gday_sun); [nlat,nlon]



  ;plot 1-map: population of selected years for map background
  yr1 = 2000
  yr2 = 2020;
  yrlen = ispan(yr1-syr,yr2-syr,1) ; 



  opt = True
  opt@PrintStat = True
  print("summary of WoE_Gday_sun after returning 99 to NA")
  cstat = stat_dispersion(WoE_Gday_sun, opt) ;detailed statistics
  print("summary of WoE_Gday_diff")
  cstat = stat_dispersion(WoE_Gday_diff, opt) ;detailed statistics
  print("summary of WoE_Gday_zero")
  cstat = stat_dispersion(WoE_Gday_zero, opt) ;detailed statistics
  ;print("summary of WoE_TWday")
  ;cstat = stat_dispersion(WoE_TWday, opt) ;detailed statistics
   


;====plot========
 imagePath = "../../image"
 imagename = imagePath+"/Fig.4"; WoE map
 wks = gsn_open_wks("pdf",imagename)
 
 plot = new(3,graphic)

 ;==set common resources
 res = True
 res@gsnDraw              = False
 res@gsnFrame             = False

 res@gsnStringFontHeightF  = 0.018;0.020
 res@tmXBLabelFontHeightF  = 0.013;0.020
 res@tmYLLabelFontHeightF  = 0.013;0.020
 res@tmYRLabelFontHeightF  = 0.013;0.020
 res@tmBorderThicknessF = 1
 res@tmXBMajorThicknessF = 1
 res@tmYLMajorThicknessF = 1
 res@tmXBMinorOn        = False           ; no lon minor tickmarks
 res@tmYLMinorOn        = False

 res@vpWidthF      = 0.8
 res@vpHeightF     = 0.3;0.4

  res@cnFillOn             = True                  ; turn on color fill
  res@cnLinesOn            = False                 ; Turn off contour lines
  res@cnLineLabelsOn       = False                 ; turn off line labels
  ;res@cnFillMode = "RasterFill"
  res@lbLabelBarOn         = True                 ;turn on/off individual label bars
  res@lbLabelAutoStride    = False ;show all labels

  ;==set res for background plot population density
  res1 = res
  res1@mpShapeMode  = "FreeAspect"
  res1@mpCenterLonF = 0.
  res1@mpMinLatF = -60.
  res1@mpMaxLatF =  75.
  ;res1@mpOutlineOn  = False; turns off the continental map outline.
  res1@mpGeophysicalLineThicknessF = 0.1    ;thiner map outlines
  res1@mpGeophysicalLineColor = "gray50"
  res1@mpOutlineBoundarySets  = "NoBoundaries"
  res1@mpOutlineSpecifiers    = "Land"
  res1@mpFillOn     = False
  res1@gsnMajorLatSpacing = 30              ; change maj lat tm spacing
  res1@gsnMajorLonSpacing = 40              ; change maj lon tm spacing

; change a default color map
  ;---Reading the default colormap into a two dimensional 254 x 4 (RGBA) array
  ;cmap = read_colormap_file("grads_rainbow")
  cmap = read_colormap_file("GMT_gray")
  cmap = cmap(::-1,:) ; reverse color map 
  ;cmap(0,:) = 1 ; set white RGBA=(1,1,1,1) for the beginning;
  res1@cnMissingValFillColor = "white"  ; missing values (default gray)
  res1@cnFillPalette        = cmap     ; use the modified color map; cnFillPalette spans the color map by default
  res1@cnFillOpacityF       = 0.6; 0.5     ; transparency of background pop

  ;-individual legend bar
  res1@lbLabelBarOn =  True;
  res1@lbOrientation          = "Vertical"
  res1@lbLabelAngleF          = 90
  res1@lbTitleString          = "Population density (per km~S~2~N~)";
  res1@lbTitlePosition        = "Right"; "Bottom"
  res1@lbTitleFontHeightF    = 0.018; 0.016
  res1@lbLabelFontHeightF    = 0.016
  res1@lbTitleDirection       =  "Across";"Down"
  res1@lbTitleAngleF          = 90
  res1@lbFillOpacityF         = res1@cnFillOpacityF
  res1@pmLabelBarWidthF       = 0.08; 0.6
  res1@pmLabelBarHeightF      = 0.3; 0.1
  res1@pmLabelBarOrthogonalPosF = 0.02; 0.10 ; move horizonal label

  res1@cnLevelSelectionMode = "ExplicitLevels"   ; set explicit contour levels
  res1@cnLevels    = (/ 1, 10, 50, 100, 500/); population density [persons/km2]

  ;--use farmer population map as the grey background
  res1@lbLabelBarOn = True
  res1@lbTitleString       = "    Population per km~S~2~N~ ~C~ doing agricultural work";
  res1@tmYLLabelsOn        = True; False
  res1@gsnLeftString        = "~F22~a"; 
  res1@gsnCenterString      = "G~S~day~N~ (dark)"; "Outdoor zero solar"; 
  plot(0)     = gsn_csm_contour_map(wks,POPR_dens,res1)

  res1@lbLabelBarOn = False;
  res1@cnFillOpacityF       = 0 ;transparent
  res1@tmYLLabelsOn        = False; True
  res1@gsnLeftString        = "~F22~b"
  res1@gsnCenterString      = "G~S~day~N~ (shade)"
  plot(1)     = gsn_csm_contour_map(wks,POPW_dens,res1); no background for c,d to clearly show WoE

  ;Outdoor sun scenario shows population density on the right
  res1@lbLabelBarOn = False;
  res1@tmYLLabelsOn        = False
  res1@gsnLeftString        = "~F22~c"
  res1@gsnCenterString      = "G~S~day~N~ (sun)"; 
  plot(2)     = gsn_csm_contour_map(wks,POPR_dens,res1)



  ;overlaying plot: show WoE as contours
  res2 = res
  ;res2@cnFillMode = "RasterFill" ; raster fill reveals more grid cells with G99>0 than default area fill
  res2@cnFillMode = "CellFill" ; cell fill shows individual grid cells more clearly than raster fill
  ;When rasterfill/cellfill is on, it is better to turn off the contour lines
  ;res@cnMonoFillPattern    = False        ; Use multiple fill patterns
  ;res@cnMonoFillColor      = True         ; Use single pattern color
  ;res@cnLinesOn            = True          ; Turn lines back on
  ;res@cnLineLabelsOn       = True;False   ; turn on line labels
  res2@cnLevelSelectionMode = "ExplicitLevels"   ; set explicit contour levels

  cmap2 = read_colormap_file("MPL_viridis"); 
  cmap2(0,:) = 1      ; set white RGBA=(1,1,1,1) for the beginning or end;
  cmap2 = cmap2(::-1,:) ; reverse color map
  ;res2@cnMissingValFillColor = "white"  ; missing values (default gray)
  res2@cnFillPalette        := cmap2; ((/1,2,3,7,8,9/),:); pick 6 colors
  res2@cnFillOpacityF        = 1; 0.8
  res2@cnNoDataLabelOn      = False; suppress the "NO CONTOUR DATA" box
 
 ;;Bottom labelbar shared by a,b 
  res2@lbOrientation        = "Horizontal"; "Vertical"    ; Rotate labelbar
  ;res2@lbLabelAngleF        = 90
  res2@lbTitlePosition       = "Bottom"; "Right"
  res2@lbTitleString         = "Warming of emergence (WoE, ~S~o~N~C)"; "Time of emergence (year)"
  res2@lbTitleFontHeightF    = 0.02; 0.016
  res2@lbLabelFontHeightF    = 0.018
  res2@pmLabelBarWidthF     = 0.7; 0.6
  res2@pmLabelBarHeightF    = 0.11; 0.1
  res2@pmLabelBarOrthogonalPosF = 0.15; 0.02
  ;res2@pmLabelBarParallelPosF = -0.08; -0.04

; If set to InteriorEdges, the labels align with the internal separators between the boxes, and there is one fewer label than the number of boxes
  ;Default option InteriorEdges means labelbar colors represent colors between levels (>= level 1 and <level 2)
  res2@cnLevels             = (/1, 1.5, 2, 2.5, 3, 3.5, 4/); , 4.5, 5, 5.5, 6/) ;WoE
  res2@lbLabelStrings       = ""+res2@cnLevels;  

  ;(default meaning of each color band is for values >= each cnLevel)
  res2@lbLabelBarOn   = False
  res2@cnInfoLabelOn  = False           ; turn off cn info label "CONTOUR FROM 1 TO 4 BY .5"
  ;overlay only works if the second plot use gsn_csm_contour, which does not have mp resources
  ;plot_0 = gsn_csm_contour(wks,WoE_TWday,res2); no valid WoE by 2099
  plot_1 = gsn_csm_contour(wks,WoE_Gday_zero,res2); zero solar; no valid WoE by 2099
  ;overlay(plot(0),plot_0); TW
  overlay(plot(0),plot_1); dark

  plot_2 = gsn_csm_contour(wks,WoE_Gday_diff,res2); diffuse;
  res2@lbLabelBarOn = True
  plot_3 = gsn_csm_contour(wks,WoE_Gday_sun,res2); full solar radiation
  overlay(plot(1),plot_2); shade
  overlay(plot(2),plot_3); sun


;=panel
  resP = True
  resP@gsnFrame         = False
  resP@gsnPanelXWhiteSpacePercent = 0; default 1 
  resP@gsnPanelYWhiteSpacePercent = 5    ;add 5% white space between plots
  ;resP@gsnPanelMainString = "Warming of Emergence (WoE) for extreme heat stress" ; set panel title
  resP@gsnPanelMainFontHeightF = 0.02
  resP@gsnPanelBottom    = 0.1          ;add some pace in the bottom
  resP@gsnBoxMargin      = 0.1          ;to accomodate the bottom legend bar

  resP@gsnMaximize         =  True ;will rotate the page
  resP@gsnPaperOrientation = "portrait"; paper orientation only works when gsnMaximize = True
 
  gsn_panel(wks,plot,(/3,1/),resP)        ; now draw as on
  frame(wks)


  print("finish WoE map") 


 ;===========plot WoE for individual models==============

 PLOT =  False
 if (PLOT) then

  ;---read WoE of each individual model, see "Data.WoE_Gday_10yravg_gridcell-ens-median-vcbc-outdoor.ncl"
  filNc  = "WoE-daytime-"+syr+fyr+"-10yravg-CMIP6-ens-12models-vcbc-outdoor.nc"
  fi = addfile(dirPath2+filNc, "r")
  WoE_Gday_sun := fi->WoE_Gday_sun
  WoE_Gday_diff := fi->WoE_Gday_diff
  WoE_Gday_zero := fi->WoE_Gday_zero
  ;WoE_TWday := fi->WoE_TWday


 imagename2 = imagePath+"/Fig.S17.Map_WoE_Gday_shade"; 
 wks2 = gsn_open_wks("pdf",imagename2)
 imagename3 = imagePath+"/Fig.S18.Map_WoE_Gday_sun"; 
 wks3 = gsn_open_wks("pdf",imagename3)
 imagename4 = imagePath+"/Fig.S19.Map_WoE_Gday_zero"; 
 wks4 = gsn_open_wks("pdf",imagename4)

  plot2 := new(12,graphic)
  plot3 := new(12,graphic)
  plot4 := new(12,graphic)
  plot5 := new(12,graphic)

  res2@vpWidthF      = 0.4
  res2@vpHeightF     = 0.15
  res2@gsnStringFontHeightF  = 0.012
  res2@tmXBLabelFontHeightF  = 0.011;0.020
  res2@tmYLLabelFontHeightF  = 0.011;0.020

  res2@mpShapeMode  = "FreeAspect"
  res2@mpCenterLonF = 0.
  res2@mpMinLatF = -60.
  res2@mpMaxLatF =  75.
  ;res2@mpOutlineOn  = False; turns off the continental map outline.
  res2@mpGeophysicalLineThicknessF = 0.5    ;thiner map outlines
  res2@mpGeophysicalLineColor = "gray50"
  res2@mpOutlineBoundarySets  = "NoBoundaries"
  res2@mpOutlineSpecifiers    = "Land"
  res2@mpFillOn     = False
  res2@gsnMajorLatSpacing = 30              ; change maj lat tm spacing
  res2@gsnMajorLonSpacing = 40              ; change maj lon tm spacing

  res2@lbLabelBarOn   = False; turn off individual labelbars
  ;shade scenario
  res2@gsnCenterString        = mo(0)
  plot2(0) = gsn_csm_contour_map(wks2,WoE_Gday_diff(0,:,:),res2);
  res2@gsnCenterString        = mo(1)
  plot2(1) = gsn_csm_contour_map(wks2,WoE_Gday_diff(1,:,:),res2);
  res2@gsnCenterString        = mo(2)
  plot2(2) = gsn_csm_contour_map(wks2,WoE_Gday_diff(2,:,:),res2);
  res2@gsnCenterString        = mo(3)
  plot2(3) = gsn_csm_contour_map(wks2,WoE_Gday_diff(3,:,:),res2);
  res2@gsnCenterString        = mo(4)
  plot2(4) = gsn_csm_contour_map(wks2,WoE_Gday_diff(4,:,:),res2);
  res2@gsnCenterString        = mo(5)
  plot2(5) = gsn_csm_contour_map(wks2,WoE_Gday_diff(5,:,:),res2);
  res2@gsnCenterString        = mo(6)
  plot2(6) = gsn_csm_contour_map(wks2,WoE_Gday_diff(6,:,:),res2);
  res2@gsnCenterString        = mo(7)
  plot2(7) = gsn_csm_contour_map(wks2,WoE_Gday_diff(7,:,:),res2);
  res2@gsnCenterString        = mo(8)
  plot2(8) = gsn_csm_contour_map(wks2,WoE_Gday_diff(8,:,:),res2);
  res2@gsnCenterString        = mo(9)
  plot2(9) = gsn_csm_contour_map(wks2,WoE_Gday_diff(9,:,:),res2);
  res2@gsnCenterString        = mo(10)
  plot2(10) = gsn_csm_contour_map(wks2,WoE_Gday_diff(10,:,:),res2);
  res2@gsnCenterString        = mo(11)
  plot2(11) = gsn_csm_contour_map(wks2,WoE_Gday_diff(11,:,:),res2);

  ;full solar
  res2@gsnCenterString        = mo(0)
  plot3(0) = gsn_csm_contour_map(wks3,WoE_Gday_sun(0,:,:),res2);
  res2@gsnCenterString        = mo(1)
  plot3(1) = gsn_csm_contour_map(wks3,WoE_Gday_sun(1,:,:),res2);
  res2@gsnCenterString        = mo(2)
  plot3(2) = gsn_csm_contour_map(wks3,WoE_Gday_sun(2,:,:),res2);
  res2@gsnCenterString        = mo(3)
  plot3(3) = gsn_csm_contour_map(wks3,WoE_Gday_sun(3,:,:),res2);
  res2@gsnCenterString        = mo(4)
  plot3(4) = gsn_csm_contour_map(wks3,WoE_Gday_sun(4,:,:),res2);
  res2@gsnCenterString        = mo(5)
  plot3(5) = gsn_csm_contour_map(wks3,WoE_Gday_sun(5,:,:),res2);
  res2@gsnCenterString        = mo(6)
  plot3(6) = gsn_csm_contour_map(wks3,WoE_Gday_sun(6,:,:),res2);
  res2@gsnCenterString        = mo(7)
  plot3(7) = gsn_csm_contour_map(wks3,WoE_Gday_sun(7,:,:),res2);
  res2@gsnCenterString        = mo(8)
  plot3(8) = gsn_csm_contour_map(wks3,WoE_Gday_sun(8,:,:),res2);
  res2@gsnCenterString        = mo(9)
  plot3(9) = gsn_csm_contour_map(wks3,WoE_Gday_sun(9,:,:),res2);
  res2@gsnCenterString        = mo(10)
  plot3(10) = gsn_csm_contour_map(wks3,WoE_Gday_sun(10,:,:),res2);
  res2@gsnCenterString        = mo(11)
  plot3(11) = gsn_csm_contour_map(wks3,WoE_Gday_sun(11,:,:),res2);

  ;zero solar (dark)
  res2@gsnCenterString        = mo(0)
  plot4(0) = gsn_csm_contour_map(wks4,WoE_Gday_zero(0,:,:),res2);
  res2@gsnCenterString        = mo(1)
  plot4(1) = gsn_csm_contour_map(wks4,WoE_Gday_zero(1,:,:),res2);
  res2@gsnCenterString        = mo(2)
  plot4(2) = gsn_csm_contour_map(wks4,WoE_Gday_zero(2,:,:),res2);
  res2@gsnCenterString        = mo(3)
  plot4(3) = gsn_csm_contour_map(wks4,WoE_Gday_zero(3,:,:),res2);
  res2@gsnCenterString        = mo(4)
  plot4(4) = gsn_csm_contour_map(wks4,WoE_Gday_zero(4,:,:),res2);
  res2@gsnCenterString        = mo(5)
  plot4(5) = gsn_csm_contour_map(wks4,WoE_Gday_zero(5,:,:),res2);
  res2@gsnCenterString        = mo(6)
  plot4(6) = gsn_csm_contour_map(wks4,WoE_Gday_zero(6,:,:),res2);
  res2@gsnCenterString        = mo(7)
  plot4(7) = gsn_csm_contour_map(wks4,WoE_Gday_zero(7,:,:),res2);
  res2@gsnCenterString        = mo(8)
  plot4(8) = gsn_csm_contour_map(wks4,WoE_Gday_zero(8,:,:),res2);
  res2@gsnCenterString        = mo(9)
  plot4(9) = gsn_csm_contour_map(wks4,WoE_Gday_zero(9,:,:),res2);
  res2@gsnCenterString        = mo(10)
  plot4(10) = gsn_csm_contour_map(wks4,WoE_Gday_zero(10,:,:),res2);
  res2@gsnCenterString        = mo(11)
  plot4(11) = gsn_csm_contour_map(wks4,WoE_Gday_zero(11,:,:),res2);

  resP2  = True
  resP2@gsnFrame         = False
  resP2@gsnPanelLabelBar = True;
  resP2@lbOrientation        = "Vertical";
  resP2@lbTitlePosition      = "Right"
  resP2@lbTitleString        = "WoE (~S~o~N~C)";"
  resP2@lbTitleDirection     =  "Across"
  resP2@lbLabelAngleF        = 0;
  resP2@lbTitleAngleF        = 90
  resP2@lbTitleFontHeightF   = 0.011
  resP2@lbLabelFontHeightF   = 0.010
  resP2@pmLabelBarWidthF     = 0.05; 0.1
  resP2@pmLabelBarHeightF    = 0.28; 0.4
  resP2@pmLabelBarOrthogonalPosF = 0.01
  ;resP2@gsnPanelRight        = 0.03 ;add some pace for the labelbar
  resP2@lbFillOpacityF        = 1

  resP2@gsnPanelXWhiteSpacePercent = 3
  resP2@gsnPanelYWhiteSpacePercent = 5    ;add 5% white space between plots
  resP2@gsnMaximize      =  True ;set this option as panel resource to maximize all plots
  gsn_panel(wks2,plot2,(/6,2/),resP2)
  gsn_panel(wks3,plot3,(/6,2/),resP2)
  gsn_panel(wks4,plot4,(/6,2/),resP2)
  frame(wks2)
  frame(wks3)
  frame(wks4)
  print("finish Fig.S17-19")

 end if



end 
